// ==================================================================================
// Shared Genomics Project MPI Codebase
// Version 1.0 30/04/2010
//
// (c) 2010 University of Manchester all rights reserved
//
// This file is distributed under the GNU General Public License, Version 2.  
// Please see the file COPYING.txt for more details
// ==================================================================================

/*!
\file
\brief	Implementation of P-value Adjustment for multiple testing. 
*/

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "adjusted.h"
#include "copyfile.h"
#include "pmaths.h"
#include "ppair.h"
#include "utils.h"
#include "suffixes.h"

struct adjusted* adjusted_create_array(int n) {
	struct adjusted* values = NULL;

	if (n > 0) {
		values = (struct adjusted*) malloc(sizeof(struct adjusted) * n);
		if (values != NULL) {
			int i;
			for (i = 0; i < n; i++) {
				values[i].bonf  =  values[i].fdr_bh = values[i].holm = 1;
				values[i].chisq = -1;
				values[i].tcnt = 0;
			}
		}
	}

	return values;
}

int adjusted_calc(struct adjusted *a, int nA, struct selected_options *ops) {	
	int i, nS = 0; // nS is elements used in the sp & schi arrays.
	double t;
	BOOL altern_pval = FALSE, jobDone = FALSE;
	struct ppair *sp = NULL;
	VECTOR_T pv_BH = NULL, pv_holm = NULL, pv_GC = NULL;

	if (a == NULL || nA <= 0 || ops == NULL) return 0;

	// See if altern_pval required, i.e. if tcnt assigned.
	// Characteristic of QT/linear algorithms.
	for (i = 0; i < nA; i++) {
		if (a[i].tcnt > 0)  {
			altern_pval = TRUE;
			break;
		}
	}

	// If have storage, do the sums.
	sp = ppair_create_array(nA);
	
	if (sp != NULL) {
		// Calculate the pvalues.
		for (i = 0; i < nA; i++) {
			if (a[i].chisq >= 0) {
				double p = altern_pval ? pT(sqrt(a[i].chisq), a[i].tcnt) : chiprobP(a[i].chisq, 1);
				if (p > -1) {
					struct ppair *pt = &sp[nS++];

					// Presume to be pvalue and original index reference to obtain.
					// original input position after data sorting.
					pt->p = p;
					pt->l = i; 
				}
			}
		}

		t = (double) nS;
		pv_BH = create_vector(nA, 1.0);
		pv_GC = create_vector(nA, 1.0);
		pv_holm = create_vector(nA, 1.0);
		
		// Works only if not zero valid tests computed.
		if (nS > 0 && pv_BH != NULL && pv_GC != NULL && pv_holm != NULL) {
			qsort(sp, nS, sizeof(struct ppair), ppair_cmp); 
			
			// Holm 
			pv_holm[0] = sp[0].p*t > 1 ? 1 : sp[0].p*t;
			for (i = 1; i < t; i++) {
				double x = (nS-i) * sp[i].p < 1 ? (nS - i) * sp[i].p : 1;
				pv_holm[i] = pv_holm[i-1] > x ? pv_holm[i-1] : x;
			}

			// BH
			pv_BH[nS-1] = sp[nS-1].p;
			for (i= nS - 2; i>=0; i--) {
				double x = (t/(double)(i+1)) * sp[i].p < 1 ? (t/(double)(i+1)) * sp[i].p : 1 ;
				pv_BH[i] = pv_BH[i+1] < x ? pv_BH[i+1] : x;
			}

			// Map the adjusted pvalues back onto the source array.
			for (i = 0; i < nS; i++) {
				struct adjusted *pv = &a[sp[i].l];
				
				pv->bonf = sp[i].p*t > 1 ? 1 : sp[i].p*t;
				pv->holm = pv_holm[i];
				pv->fdr_bh = pv_BH[i];
			}

			jobDone = TRUE;
		}
	}
	if (sp != NULL) free(sp);
	if (pv_BH != NULL) free(pv_BH);
	if (pv_GC != NULL) free(pv_GC);
	if (pv_holm != NULL) free(pv_holm);
	
	if (jobDone) return 1;
	else return 0;
}

int adjusted_write(struct adjusted *a, int nA, struct selected_options *ops) {
	BOOL jobDone = FALSE;
	FILE *out = NULL;

	if (a == NULL || nA <= 0 || ops == NULL) return 0;
	else if (strlen(ops->szAdjustedFileName) == 0) return 0;

	out = fopen(ops->szAdjustedFileName, "wb");
	if (out != NULL) {
		if (fwrite(&nA, sizeof(int), 1, out) == 1) {
			if(fwrite(a, sizeof(struct adjusted), nA, out) == nA) jobDone = TRUE;
		}
		fclose(out);
	}

	if (jobDone) return 1;
	else return 0;
}

struct adjusted* adjusted_load(struct selected_options *ops, int *size) {
	int _size = 0;
	FILE *in = NULL;
	struct adjusted *data = NULL;

	if (ops == NULL || size == NULL) return data;
	else if (strlen(ops->szAdjustedFileName) == 0) return data;

	*size = 0;
	in = fopen(ops->szAdjustedFileName, "rb");
	if (in != NULL) {
		if (fread(&_size, sizeof(int), 1, in) == 1) {
			if (_size > 0) {
				data = (struct adjusted*) malloc(sizeof(struct adjusted) * _size);
				if (data != NULL) {
					if (fread(data, sizeof(struct adjusted), _size, in) == _size) {
						*size = _size;
					} else {
						free(data);
					}
				}
			}
		}
		fclose(in);
	}

	return data;
}

int adjusted_write_and_copy(struct adjusted *a, int nA, struct selected_options *ops, int rank) {
	char remoteFilepath[FILE_PATH_LENGTH];
	BOOL jobDone = FALSE;
	char tmp_szRunId[TOKEN_LENGTH];

	if (a == NULL || nA <= 0 || ops == NULL) return 0;

	strcpy(tmp_szRunId, ops->szRunId);
	sprintf(ops->szRunId, "%s_%i", tmp_szRunId, rank);

	sprintf(ops->szAdjustedFileName, "%s%s%s", ops->ldir, ops->szRunId, SUFFIX_PADJ_FILE, rank);
	sprintf(remoteFilepath, "%s%s%s", ops->rdir, tmp_szRunId, SUFFIX_PADJ_FILE);

	if (adjusted_calc(a, nA, ops)) {
		if (adjusted_write(a, nA, ops)) if (copyAndDeleteFile(ops->szAdjustedFileName, remoteFilepath)) jobDone = TRUE;
	}

	strcpy(ops->szRunId, tmp_szRunId); // Restore the original Run ID

	if (jobDone) return 1;
	else return 0;
}

struct adjusted* _adjusted_load(struct selected_options *ops, int *size) {
	if (ops == NULL || size == NULL) return NULL;

	sprintf(ops->szAdjustedFileName, "%s%s%s", ops->rdir, ops->szRunId, SUFFIX_PADJ_FILE);
	return adjusted_load(ops, size);
}